function z = h(x,y,k,nu)
switch k
case 1  % for stokes problem
    z = 18*(y.^2 - x.^2);
case 2  % for stokes problem
    z = -48*nu + 96*(y.^3.*x-x.^3.*y);
case 3  % for stokes problem
    z = -4*nu*sin(x+y);
case {0,4,5,105,7,107,10,11,15,18,33,66,88,158,200,150,300,301,104}
    z = zeros(size(x));
case 6
    z = -12*y.*(x-1).^3;
case 8
    z = 288*x.^2 - 384*x + 192 -12*y.^4 +48*y.^3 -192*y - 96*x.^3 +12*x.^4 - 48*nu +24*x.^3 .*y ...
        -144*x.^2 .*y + 288*x.*y - 24*x.*y.^3;
case 9
    z = x.*y -2*x.*y.^2 + 4/3*x.*y.^3;
case 12
    z = (1-x).*(x.^2.*(12*exp(2*y)+nu*exp(y))-x.*(24*exp(2*y)+2*nu*exp(y))+12*exp(2*y)+13*nu*exp(y));
case 13
    z = -3/2*nu*(sin(y+x-1)+sin(-y+x-1))+25/2*nu*(sin(3*x-3+y)+sin(3*x-3-y))+3/4*(cos(-2*y+2*x-2)-...
        cos(2*y+2*x-2))+3/8*(cos(4*x-4+2*y)-cos(4*x-4-2*y));
case {101,106,109}
    z = zeros(size(x));
case 400
    z = -48*nu*ones(size(x));
case 500
    z = -1680*nu*(x.^4+y.^4);
end;
